The role of symmetries in adiabatic quantum algorithms 
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Exploiting the similarity between adiabatic quantum algorithms and quantum phase transitions, 
we argue that second-order transitions - typically associated with broken or restored symmetries 
- should be advantageous in comparison to first-order transitions. Guided by simple examples we 
construct an alternative adiabatic algorithm for the NP-complete problem Exact Cover 3. We show 
numerically that its average performance (for the considered cases up to 0{2O} qubits) is better 
than that of the conventional scheme. The run-time of adiabatic algorithms is not just determined 
by the minimum value of the fundamental energy gap (between the ground state and the exited 
states), but also by its curvature at the critical point. The proposed symmetry-restoring adiabatic 
quantum algorithm only contains contributions linear and quadratic in the Pauli matrices and can 
be generalized to other problem Hamiltonians which are decomposed of terms involving one and 
two qubits. We show how the factoring problem can be cast into such a quadratic form. These 
findings suggest that adiabatic quantum algorithms can solve a large class of NP problems much 
faster than the Grover search routine (which corresponds to a first-order transition and yields a 
quadratic enhancement only). 

PACS numbers: 03.67.-a, 03.67.Lx, 73.43.Nq. 



I. INTRODUCTION 



A. Adiabatic quantum algorithms 



The overwhelming potential of quantum computers in 
number factorization [l| and database search [2| has ini- 
tiated a lot of research. In the conventional (sequen- 
tial) quantum algorithmic approach, unitary single and 
two-qubit operations are applied on a quantum state 
and the result of the computation can then be extracted 
(probabilistically) by measurement. Unfortunately, the 
promises of quantum computing are strongly hampered 
by the devastating influence of decoherence: Inevitable 
couplings with the reservoir tend to destroy the frag- 
ile quantum coherence. The available error correction 
schemes require further auxiliary qubits [3| and - al- 
though even strongly correlated errors can in principle 
be corrected 0] _ become increasingly complex with ad- 
mitting more and more general errors. 

Farhi et al. have proposed an alternative scheme with 
an inherent error protection: adiabatic quantum compu- 
tation [1, @. In this scheme, the solution to a problem 
is encoded as the (unknown) ground state of a (known) 
problem Hamiltonian Hp, which is separated by a suf- 
ficiently large energy gap from excited states. Compu- 
tation will start with a different initial Hamiltonian Hi 
possessing a known and easily preparable ground state, 
which is also well separated energetically from all other 
states. The system is prepared in this ground state and 
then the initial Hamiltonian Hi is continuously deformed 
into the final problem Hamiltonian Hp. The adiabatic 
theorem guarantees that if the evolution is slow enough, 
the system will stay near its instantaneous eigenstate 
and will thus finally reach the ground state of the prob- 



lem Hamiltonian Q. Thus, the solution to the problem 
can be extracted by measurement of the final quantum 
state. A nice advantage of this quantum computation 
scheme (which is believed to be polynomially equiva- 
lent to sequential quantum computation Q) lies in the 
robustness of the ground state against the influences 
of decoherence - a sufficiently cold reservoir provided 
0, M EH HI IH M EE M M, E3 : The ground state 
cannot decay and phase errors do not play any role, i.e., 
errors can only result from excitations. The speed (total 
evolution time T) at which the interpolation between the 
two Hamiltonians can be performed without strongly ex- 
citing higher eigenstates of the system is related to the 
inverse gap between the instantaneous lowest and first 
excited eigenvalues @. The total evolution time T re- 
quired to reach a fixed fidelity in the final state can be 
interpreted as the algorithmic complexity of the adiabatic 
computation scheme. With using a constant deformation 
speed (ignoring the structure of the fundamental energy 
gap) the adiabatic run-time will scale as T = 0{g^? n }, 
whereas with knowledge on the fundamental gap this can 
under fairly moderate assumptions even be improved to 
T = 0{g~l} 00. 

This also highlights the main obstacle in adiabatic 
quantum computation: Typically, the minimum funda- 
mental gap g m i n decreases strongly with increasing sys- 
tem size (number of qubits) for nontrivial problems. For 
example, in the adiabatic version of Grover's algorithm 
[2l| . the minimum fundamental gap decreases exponen- 
tially <7 m i n = 2~ n / 2 with system size n. An analogous 
exponential scaling of the minimum fundamental gap has 
been found in other adiabatic quantum algorithms with 
similar initial Hamiltonians |22|]. However, it has been 
argued that this exponential scaling is a result of the un- 
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favorable choice of the initial Hamiltonian and can be 
avoided for more suitable choices [23|]. In these cases, 
however, Omi n is not known analytically. Numerical anal- 
ysis in 0, |(| seems to indicate a quadratic scaling T ~ n 2 
of the algorithmic complexity with the system size n - but 
this favorable scaling could perhaps be just a consequence 
of choosing particularly simple problems, cf. [2J, [25[. 
Therefore, the speed-up attainable with these adiabatic 
quantum algorithms is still an open question. 



B. Quantum phase transitions 

As has been noted earlier 0, [2tJ , adiabatic quantum 
algorithms display a remarkable similarity with sweeps 
through quantum phase transitions: During the adia- 
batic interpolation, the ground state changes from the 
simple initial ground state of Hi to the unknown solu- 
tion of some problem encoded in Hp. Typically, on the 
way from H\ to Hp, one encounters a critical point where 
the fundamental gap (which is sufficiently large initially 
and finally) becomes very small. Near the position of the 
minimum gap, the ground state will change more drasti- 
cally than during other time intervals of the interpolation 
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n>0 



E n — E 



(1) 



where we have inserted an identity and expressed 
^V'olV'n^ by the time-derivative of the eigenvalue 

equation H(t)\ip (t)) = E (t) \ipo(t)) . If the rele- 
vant fundamental energy gap (the smallest one with 
(ipo\ H \ip n ) ^ 0) scales inversely with the system size, we 
see directly that in the infinite-size limit, the ground state 
will change non- analytically at the critical point. This 
singularity would completely prohibit an adiabatic evolu- 
tion in this limit. However, for practical problems one is 
interested in the finite-size scaling of an algorithm, where 
it makes a huge difference whether its computational 
complexity increases exponentially or merely polynomi- 
ally with the system size n. The same principle applies to 
the efficient suppression of thermal excitations when the 
quantum system is coupled to a low-temperature bath 
[9(, where it may prove experimentally difficult to apply 
polynomially small temperatures and most likely infeasi- 
ble to apply exponentially small temperatures. 

In the following, we shall exploit the analogy between 
adiabatic quantum algorithms and quantum phase tran- 
sitions further in order to gain additional insight into 
these issues. Beyond the behaviour of the minimum gap, 
the complexity of implementing the involved Hamiltoni- 
ans (number of interactions between different qubits) is 
of experimental importance and shall also be addressed. 



II. MOTIVATION 



Grover Search Routine 



The adiabatic version of Grover's search al- 
gorithm is defined by the linear interpolation 
H(s) = (1 — s) Hi + sHp, where s e [0,1] between 
the Hamiltonians I2l| 



H l = 1 - \S) (S\ , 
where the initial ground state 



H, 



N-l 



\S) 



(2) 



(3) 



z=0 



N 



is the superposition vector of all 
states in the computational (o~ z ) 

HH(|o) + |i»A/2 = (U) + |T))A/2. 

the final ground state \w) is some 
computational basis state. 

It is straightforward to derive the instantaneous spec- 
trum of H(s) by Erhard-Schmidt orthogonalization, for 
example. One obtains the fundamental gap 
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basis and 
In contrast, 
distinguished 



g(s) = 



1 i i - ) - 



(4) 



■o/i 00 = \ 

1. Thus, the adiabatic 



[l±g(s)} as 



and the two non-trivial levels E, 
well as E 2 (s) = . . . = E N -i(s) 
Grover algorithm continuously transforms the initial vac- 
uum \S) towards the final ground state \w), see also fig- 
ure [TJ From Equation ([¥]) we may infer an exponentially 
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FIG. 1: Sketch of the two lowest energy levels for the Grover 
problem as a prototypical example for a first-order quantum 
phase transition in the infinite size limit n | oo (where the 
energy gap vanishes and the avoided level crossing becomes a 
real crossing). 



small minimum fundamental gap g mm = yl/2" at the 
critical point s cr ;t = 1/2. 

In the infinite-size limit (n — > oo), the first derivative 
dEo/ds of the ground state energy E {s) jumps at the 
point of the phase transition s cr ; t = 1/2. In view of the 
equality 



dE (s) 
ds 



V>o(s) 



dH{s) 



ds 



■00 (s) 



(5) 
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this discontinuity goes along with a jump in certain ex- 
pectation values (e.g., (dH/ds)) and hence order param- 
eters such as the horizontal magnetization (a x ) change 
abruptly at the critical point s cr j t = 1/2. Conventionally, 
this property (discontinuous order parameter) is used to 
classify the phase transition as a first-order transition. 

From Eqn. ([I]) it follows that the ground state basically 
jumps from |-0 o (s < 1/2)) = \S) to \ip (s > 1/2)) = \w). 
Such an abrupt change is a general feature of first-order 
quantum phase transitions and can be understood in 
terms of a time-dependent energy landscape. Assum- 
ing a separable state (compare also [28|, |29( for a similar 
approach) 

n 

k) = 0[cos(^)|O) +sin(^) |1>], (6) 

8=1 

and re-ordering the states in the computational basis 
such that \w) = |1 . . . 1) (the initial state ([3]) is invariant 
to this transformation) we can calculate the semiclassical 
energy landscape 

E(s,<p) ee (<p\H( 8 )\<p) 

_ [cos(y) +sin(y)f" 

' 2™ 
SBm n (<p), (7) 

shown in Fig. [2j We see that even in this semiclassical 



The time-dependent Hamiltonian ([2} does not have a 
conserved quantity, such that the initial ground state 
| S) does not break any symmetry of the initial Hamil- 
tonian. We conjecture that in the absence of additional 
conserved quantities (related to continuous symmetries of 
the Hamiltonian), there are always two (or more) com- 
peting local minima whose energies depend on the inter- 
polation parameter s. At the critical point s cr it, the or- 
der of the s-dependent energy minima changes and thus 
the global minimum makes a jump. Consequently, in 
order to stay in the true ground state, the quantum sys- 
tem has to tunnel through the barrier depicted in figure 
[2j Since one would naturally expect that all the energy 
scales and hence also the size of the tunneling barrier in- 
crease roughly linearly with system size n, the tunneling 
probability decreases exponentially with the number of 
qubits n. This intuitive picture suggests that the run- 
time T needed to stay in the ground state scales expo- 
nentially for all first-order transitions. Indeed, since the 
tunneling amplitude provides the coupling between the 
competing ground states and hence determines the min- 
imum gap 5mi n in their avoided level crossing, such an 
exponential scaling is precisely what one finds for the adi- 
abatic Grover algorithm and further adiabatic algorithms 
[22l I23I l2l| with similar initial Hamiltonian - which all 
correspond to first-order transitions. 



B. Transverse Ising model 




FIG. 2: Semiclassical energy landscapes (0 for the Grover 
problem with n = 10 spins for different interpolation param- 
eters a. In order to stay in the global minimum, the sys- 
tem has to tunnel through an energy barrier that does not 
vanish throughout the interpolation, i.e., the time-dependent 
position (red dotted lines) of the two local minima is not 
connected. Circle symbols show the position of the global 
minimum. 



After having discussed first-order transitions, let us 
turn to a prototypical model for a second-order quan- 
tum phase transition [301 ] . The one-dimensional trans- 
verse Ising model is given by linear interpolation between 

n n 

Hi = -yv* , h f = -^2<jf(ji +ll (8) 

e=i 1=1 

where u x J v denote the Pauli spin matrices acting on the 
£ th qubit and periodic boundary conditions of i+1 = o\ 
are assumed. At any point during the interpolation, 
the Hamiltonian H(s) can be diagonalized analytically: 
The successive application of Jordan- Wigner, Fourier, 
and Bogoliubov transformations [30] map the interact- 
ing spin-l/2-Hamiltonian |5]) to a set of non-interacting 
fermionic quasi-particles 

B = fa* ~\) > W 

with fermionic creation and annihilation operators 7I, 7fc 
and single quasi-particle energies 



picture the two minima are always separated - either by a 
potential barrier or by large energy differences. This cor- 
responds to the existence of a tunneling barrier between 
the two minima. 



e k (s) = 2y/l - 4cos 2 (fca/2)s(l - s) , (10) 

see figure [3] The wavenumber ka covers the range 
ka = (1 + 2Z) ir/n : \ka\ < n. Both Hamiltonians in 
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FIG. 3: Sketch of the two lowest energy eigenvalues for the 
Ising model as a prototypical example of a second-order quan- 
tum phase transition. In contrast to figure [T] the geometry 
near the minimum gap does not correspond to an isolated 
avoided crossing. 



([8]) obey a 180-degree rotational symmetry around the 
of/-axes (bit-flip) 



n 




n 






Hp,® of 






e—i 



(11) 



which transforms of to —of. Since due to Eqn. (TTTj) 
the even bit-flip parity of the unique initial ground state 
\S) is a conserved quantity, the symmetry of the two- 
fold degenerate final ground state |\Pq) = || • • • t) an d 
|*o) = II • ■ ■ I) i s broken, such that for adiabatic evolu- 
tion, the system will end up in the macroscopic superpo- 
sitions (Schrodingers cat) state 



l*o von > 



T) + U ■■■!) 

V2 



(12) 



The minimum gap can be obtained from equation (flT 
and scales polynomially g m i n = 0{l/n}. Therefore, the 
adiabatic runtime T does also scale polynomially - a 
constant speed interpolation, for example, yields T ac i = 
0{n 2 } [31(. Furthermore, since the wavenumber ka 
covers the range ka = (1 + 2Z) n/n : \ka\ < n, the 
infinite-size limit of the ground state energy can be ob- 
tained by replacing the sum over the single-particle en- 
ergies by an integral 

E^is) = - lim 



-^(2v^(l- 



«)) , (13) 



where £{x) denotes the complete elliptic integral. In con- 
trast to the previous subsection III Al the energy density 
£(1) = 1 and its first derivative £'{1) = are well defined 
at the critical point. Since (H(s)) and (dH/ds) are con- 
tinuous across the transition, there is no jump in order 
parameters such as (a x ) and the ground state changes less 
abruptly. However, the second derivative of the ground 
state energy diverges at the critical point £"(l) — oo and 
hence, we classify the Ising model as a second order quan- 
tum phase transition, compare also [30(. Consequently, 
we have a symmetry-breaking quantum phase transition 
of second order. 



As in the previous section, we can derive a semiclassi- 
cal time-dependent energy landscape with the separable 
ansatz ^ to obtain 



E, 



,(s,ip) = nj - (1 - s) sin(2y>) 

- S [cos 2 M-sin%)] 2 }, (14) 

where the bit-flip invariance (jTTJ) is reflected by a mirror 
symmetry of the energy landscape around tp — tt/4, see 
figured] We see that symmetry-breaking/restoring quan- 
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FIG. 4: Semiclassical energy density landscape for the Ising 
model in a transverse field for different interpolation parame- 
ters s. In contrast to first-order transitions, there is no energy 
barrier between the initial and the final ground state, i.e., the 
time-dependent positions (red dotted line) of the minima are 
always connected. 

turn phase transitions are typically of second (or higher) 
order, i.e., the ground state does not change as abruptly 
as in first-order transitions. Consequently, the system 
does not need to tunnel through a barrier in order to 
stay in the global vacuum. As a result, the reason for 
the exponential scaling of the run-time T we found for 
first-order transitions is absent in the case of quantum 
phase transitions of second or higher order, which is con- 
sistent with thepolynomial scaling of gap [3(| and adia- 
batic runtime [31] . 

Evidently, the transverse Ising model can also be con- 
sidered an adiabatic quantum algorithm, although of 
course the final ground states arc trivial the problem 
to be solved is simply "How can all bits have the same 
value?" Note that, reversing the evolution and slowly 
interpolating from H-p towards H\ in ([5]), the final state 
will be close to the true unique ground state IS*) if and 
only if initialized with the superposition (fT2"]l . 



C. Mixed Case 

Even though the concept of broken or restored sym- 
metries is very useful for classifying quantum phase 
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transitions, it should be stressed that symmetry break- 
ing/restoration alone does not guarantee a second-order 
transition. As an intuitive counterexample, one may con- 
sider the linear interpolation between 



\S) (S\ , 



n 1 

E^ 1 



(15) 



which has an equivalent (up to shifting and scaling to 
obtain positive definiteness) final Hamiltonian as |(5J| but 
differs in the initial Hamiltonian, for which we have cho- 
sen the initial Hamiltonian of the Grover problem [[2]). 
Evidently, the ground state symmetry of the final Hamil- 
tonian is also broken, since the bit-flip parity is also con- 
served. However, an analysis of the semiclassical energy 
landscapes generated with a separable ansatz ([6]) yields 



E n (s,<p) = (l-s)U 



[cos(yj) + sin(y»)] 



In 



+ y {l-[cos 2 M-sin 2 M] 2 }, (16) 




0.2 0.4 0.6 

interpolation parameter s 

FIG. 6: Lowest eigenvalues of the Hamiltonian (|15[) for 10 
(bold black lines) and 15 (thin red lines, shifted for better 
readability) qubits in the relevant subspace which is even un- 
der bit-flip. One can clearly see that the spectrum displays an 
avoided level crossing at the critical point - thus correspond- 
ing to a first-order transition. 



which also exhibits a mirror symmetry at (p — 7r/4 but 
always has a tunneling barrier between the vacua, see 
figure [5] The level structure displays the geometry of an 
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FIG. 5: Semiclassical energy landscapes for a symmetry- 
breaking quantum phase transition - which is, however, not of 
second but of first order. In spite of the symmetry breaking, 
there is a tunneling barrier throughout the interpolation and 
a jump between the initial and the final ground state(s), i.e., 
the position of the vacua (dotted red lines) is not connected. 

avoided level crossing at the critical point, i.e., it corre- 
sponds to a first-order phase transition, see Figure O In 
accordance with our previous arguments, it can be shown 
analytically that the gap scales exponentially in such a 
situation ^HHH. 

In order to understand why the phase transition is of 
first order, it might be useful to recall that the initial 
Hamiltonian is a projector that involves n-qubit interac- 



tions 

fl I = l-|S)(S| = l-0l±fl. (17) 

The bit structures of the initial and the final Hamiltonian 
are very different and, in this sense, their "distance" in 
the space of all matrices is very large. 

D. Conjecture 

We have seen for the Grover example that first order 
quantum phase transitions are associated with an expo- 
nential scaling of the adiabatic runtime with the system 
size n. For local Hamiltonians, this observation can be 
generalized [33j]: For a typical first order transition one 
has two locally distinguishable competing ground states 
\w) and \s) that exchange their energetic favorability 
right at the critical point (w\ H \w)\ a = (s\ H \s)\ s . . 
Their overlap is exponentially small (due to local dis- 
tinguishability) and for a local reasonable (with polyno- 
mial matrix elements) Hamiltonian this also implies that 
(s\H\w) cx exp{— 0(n)} is also exponentially small. In 
the two-dimensional subspace spanned by these states the 
energy gap equates to g cr i t = 2\(s \ H \ w) |, which becomes 
therefore exponentially small. 

Our conjecture is therefore that for adiabatic quantum 
algorithms 

a. ) higher order quantum phase transitions are more ad- 
vantageous than first order quantum phase transitions 
and 

b. ) such higher order transitions might be induced by 
interpolations respecting a conserved quantity (leading 
to spontaneous symmetry breaking in the degenerate 
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ground state of one Hamiltonian) in combination with 
similar bit structures in initial and final Hamiltonians. 

Of course, these considerations are very intuitive and 
by no means conclusive - but they will hopefully help us 
to design more powerful quantum algorithms. For exam- 
ple, it would of course be nice to have no phase transition 
at all (which corresponds to a constant lower bound on 
the fundamental energy gap [32j]) but numerous numer- 
ical [H, [H, [13, H [35| and analytical [H, H| evidence 
suggests that this is not the case for hard optimization 
problems. 



III. RUNTIME SCALING FOR LINEAR 
QUENCHING 

The main measure for the computational complexity of 
adiabatic quantum algorithms is the minimum runtime 
T. In the following, we derive an estimate of the runtime 
T based on the level structure of the adiabatic quantum 
algorithm. We will assume a linear quench s = t/T and 
a single isolated position of the minimum energy gap. In 
cases where the position and size of this minimum energy 
gap is explicitly known, one can improve the runtime by 
adapting the interpolation speed (i.e., moving fast away 
from the critical point and slow in its vicinity) p^ . l2lf . 
However, since the exact position of the minimum gap 
is not known a priori for most of the interesting cases - 
although its approximate position can be estimated by 
perturbative and other (such as e.g., Gershgorin's circle 
theorem) methods Q - we shall work with a constant- 
speed interpolation. 

As derived in [l9| . the final occupation amplitude of 
the first excited state ai(s = 1) can be obtained formally 
from the adiabatic expansion via 



9(s) 



x exp I -iT j g(s')ds' 



(18) 



where Fq\ (s) denotes the transition matrix element 
of H'(s) in the instantaneous energy basis, ao(s) the 
ground-state amplitude, and -fi is a pure phase (includ- 
ing the Berry phase). The ds- integration along the real 
axis in the above integral can be deformed in the lower 
complex half- plane [ljj], where it is visible that for slow 
interpolations (large T), the outer integrand in (fP8|) is ex- 
ponentially suppressed (adiabatic approximation). How- 
ever, such a deformation will be limited by singularities 
of g (s) at s located in the complex plane near the min- 
imum of g(s) on the real axis. This determines how the 
adiabatic runtime T has to scale with the spectral prop- 
erties in order to suppress the excitation amplitude in 



IS) efficiently 



ft 



3?(s)+i3(s)/2 

iT J g{s)ds 
b 



»1, 



(19) 



for a detailed discussion see [T{| . Assuming a sufficiently 
smooth behavior, we Taylor-Laurent expand the funda- 
mental gap near its minimum (which is a saddle-point in 
the complex plane) 



9(s) 



~r~ C m j n (s ^crit) j 



(20) 



where g m i n denotes the value and c m i n the curvature of 
the fundamental gap at the critical point s = s cr i t . The 
singularities of \/g(s) are then found at approximately 
s s cr it ± i\J <7min/cmin, which yields for the left hand 

side of era 



ft 



Scrit 3: 



a/2 



iT 



g(s)ds 



OlT 



3/2 
,1/2 



(21) 



Consequently, the runtime necessary to suppress the ex- 
citation amplitude (ThS) efficiently scales as 



T = 



g 3 - 

■J mm 



(22) 



Evidently, the runtime does not only depend on the value 
of the minimum gap g m i n but also on its curvature c m i n 
at the critical point s = s cr it- At first sight, this may 
seem counterintuitive, since from a naive interpretation 
of the adiabatic theorem (focusing on the minimum gap 
only), one would expect that a small curvature (which 
implies a longer persistence of a small gap) should lead 
to longer run-times. However, it should be kept in mind 
that a large curvature means a rapid change of the spec- 
tral characteristics of the Hamiltonian and hence the sys- 
tem will find it harder to follow the evolution in order to 
stay in the ground state [l^, [2(| • 

As a consistency check, we show that (|2"2"|) correctly 
reproduces the scaling of the adiabatic runtime found 
earlier for the Grover problem and the Ising model in case 
of constant-speed interpolation: For the Grover model, 
one obtains from equation ([?) 



T 



'5c(£crit) 



2N-1 + 



(23) 



which implies an exponential scaling of the adiabatic run- 
time (since N = 2 n ) with the system size, compare also 
fTg| . |2"H and subsection III Al Likewise, one obtains for 
the Ising model from equation (jTDJ) 



T 



/ jjfCgcrit) 
5l 3 (Scrit) 



1 



1 



(24) 
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a quadratic scaling of the adiabatic runtime, see also [3l| 
and subsection III Bl 

We have numerically calculated the curvature and 
value of the minimum gap for interesting adiabatic op- 
timization problems (see the following subsections) and 
have then related the adiabatic runtime T with its esti- 
mator v/ Cmin/fl'miii - Our numerical experiments also con- 
firm this scaling law quite nicely, see figure [7) 
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FIG. 7: [Color Online] Scaling of the adiabatic runtime from 
figure [8] versus a/ Cmin/fl^in ^ or different adiabatic algorithms. 
All data sets agree nicely with the relation (dashed fit lines) 
predicted by equation (|22[1 . With the exception of the Grover 
algorithm, each symbol corresponds to the median of 100 
random instances with a unique satisfying agreement with- 
out constraints (brown) on the number of clauses or with 
m < 2/371 clauses (black and red). Values for the linear fits 
(dashed lines) are given in brackets. Horizontal error bars 
give the 99% confidence interval on the median, whereas ver- 
tical error bars result from the finite precision when deter- 
mining the necessary runtime by repeated integration of the 
Schrodinger equation. The data sets for Grover interpola- 
tion stem from fl9l ]. where a different success probability of 
.Pflnai = 0.75 had been demanded. 



IV. QUADRATIC SPIN HAMILTONIANS 

Our conjecture (subsection III Dp can be exploited to 
improve the search for ground states of quite general 
quadratic Hamiltonians. ft is well known that the ground 
state of a frustrated anti-ferromagnet under the influence 
of a local field as described by the Hamiltonian 

n n 

flp = 5>*f+ £ ho\a) (25) 

z=l i,j = l 

is in general extremely (NP-)hard to hnd [46]. In above 
equation, the coefficients fi and /y = fa (with fa = 0) 
denote the local held and the interaction topology, re- 
spectively. More technically speaking, deciding whether 
the minimum energy of (|25p lies below a certain thresh- 
old belongs to the class NP, i.e., it can be verified for 



a given solution with polynomial effort - whereas actu- 
ally finding the solution can be much (e.g., exponentially) 
harder. For Exact Cover 3, this decision problem is even 
NP-complete @ , which implies that all other problems in 
NP (i.e., decision versions of factoring, traveling salesmen 
problem, etc.) can be mapped to it with using only poly- 
nomial overhead. There is no classical algorithm known 
that delivers a solution using only polynomial resources 
in n and the existence of such an algorithm would dis- 
prove the commonly accepted conjecture P ^ NP, which 
is of vital importance in classical complexity theory. 

The discovery of a polynomial-time algorithm would 
have enormous consequences: The (worst-case) perfor- 
mance of all known classical algorithms scales exponen- 
tially in the problem size n and already the solution of 
relatively small problems becomes classically impossible 
for all practical purposes. This is for example exploited 
in encryption protocols, see also section [VI] and section 
IVIIIl in the appendix. 

A. XVZ-Network 

In order to find a good adiabatic quantum algorithm to 
reach the ground state of (|2"5|) , we want to use the insight 
from the previous sections and try to construct an initial 
Hamiltonian H\ such that the interpolation from Hi to 
Hp corresponds to a phase transition of second or higher 
order. Since broken or restored symmetries are typically 
(though not always) associated with second-order phase 
transitions, we shall demand that Hi shares a symmetry 
with Hp and that this symmetry is broken in the initial 
ground state but restored in the final solution state. Fur- 
thermore, Hi should retain the interaction topology (i.e., 
bit structure) of Hp. 

To this end, one can exploit the evident symmetry of 
the problem Hamiltonian (|25p with respect to rotations 
around the axis generated by 

n 

£* = ^>f, (26) 
»=i 

which we will further-on denote as Hamming-weight op- 
erator, since its eigenvalues are associated with the Ham- 
ming weight of the corresponding subspace. An obvious 
example for a Hamiltonian where this rotational symme- 
try is spontaneously broken is the SO (3)-invariant ferro- 
magnetic Hamiltonian 

n 

11 \" = - J2 l/yl°"i ,<T J 

n 

= - E \h\[*f^ + + *?^] , (27) 

where we have used the same interaction topology fa as 
in (|25|) . Clearly, this Hamiltonian is invariant under ro- 
tations around arbitrary axes - including (|26p - whereas 
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its (n + l)-fold degenerate ground state singles out one 
specific direction, for example ||t • • ■ T)< 

This degeneracy immediately poses the question of 
which of these ground states should be used as the ini- 
tial state. In analogy to the Ising model, where only 
the subspace with even bit-flip parity was relevant, this 
question can be answered by the conserved quantity (|26[) . 
Since the final solution state will have a fixed Hamming 
weight A, we should start in the corresponding subspace 

Ha = {|V>) : S z = A |^)}. Obviously, |TT T) an d 

1 I J. ... I) correspond to A = n and A = — n, respectively, 
and are not suitable (except in trivial cases). However, 
the state \S) = |— >— > •••—>) contains a superposition of 
all Hamming weights and can be used to project out an 
initial state with any desired value of Afe = 2k — n via 



iiifc) 



\ 



2n 



2n 
k=0 



exp 



2m ■ 



A* 



2n+ 1 



k |S>(28) 



Up to normalization, above formula is just the Fourier 
decomposition of the Kronecker symbol <5(A — E z ) and 
involves single-qubit rotations only. The state (|28[) can 
be prepared efficiently by different approaches such as 
for example projective measurements or adiabatic evolu- 
tion [3(| . Alternatively, one could use relaxation with an 
appropriate energy penalty in the Hamiltonian such as 
(E z -A) 2 . 

Of course, in order to prepare the initial state |in) cor- 
rectly, one has to know the Hamming weight A of the 
solution bit-string - and in most adiabatic algorithms, 
this knowledge will not be available in advance. But 
since there are for an n-bit problem only n + 1 subspaces 
with different Hamming weights, a naive testing of all 
these possibilities corresponds to a polynomial overhead 
only. Moreover, for many problems such as exact cover 
3 discussed in the next section, one can guess the rough 
value of A and thereby limit the number of trials. 



B. XY-Network 

So far, we considered the fully SO (3) rotationally sym- 
metric ferromagnetic Hamiltonian (|27p . However, nu- 
merically we have found (at least for the exact cover 3 
problem) that the planar ferromagnetic model 



(29) 



with the same interaction topology /y yields a better 
performance of the adiabatic algorithm (in average). In 
contrast to the spherically symmetric Hamiltonian (|27p, 
which is invariant under rotations around an arbitrary 
axis, this planar Hamiltonian is merely axially symmet- 
ric, i.e., invariant under rotations around the E z -axis, 
i.e., = 0. 



Our conjecture is that for cases where the H* yz initial 
Hamiltonian performs significantly worse, not only E z 
is conserved exactly but - since this initial Hamiltonian 
commutes with n ■ S = n er^ - one might obtain ad- 
ditional nearly conserved quantities via [n ■ S, Hp] rs 0. 
Since for an arbitrary vector n, the final ground state 
may not have the symmetry, such nearly conserved quan- 
tities would strongly hamper the evolution from the ini- 
tial ground state to the final one. Of course, the prob- 
ability to find such an axis n is greatly reduced for the 
Hj V initial Hamiltonians, which is consistent with our 
observation of increased performance. 

Furthermore, the state \S) = |— ►— > • • • — ►) (as well as 
arbitrary rotations of it around the E z -axis) is no longer 
the exact ground state of (f2"9")l in general. These states 
are separable (no entanglement) whereas the true ground 
state of (|29| is already entangled and cannot be given 
analytically for the general case. Fortunately, the mean- 
field approximation typically works reasonably well for a 
large number n ^> 1 of pseudo-randomly connected spins 
and hence the state \S) should provide a good approx- 
imation to the ground state. Indeed, numerically [13] 
we found that the state in (f2"B"|) possesses a large overlap 
(over 90%) with the exact ground state (in the relevant 
subspace Ha) and therefore provides a good initial state 
for the computation. Thus the exact SO (3) degeneracy 
of the ground state in the rryz-network (|27[) is replaced by 
an approximate 0(2) degeneracy (mean-field approxima- 
tion) in the zy-network (|2"9"|) . Nevertheless, both models 
imply a symmetry-restoring transition and thus should 
be of second order. Note that interestingly, neither |27p 
nor (|29p are diagonal in the Hadamard basis, such that 
the arguments in [25| do not directly apply. In addition, 
these Hamiltonians differ strongly from projection oper- 
ators (compare the discussion in [22], HH ) . 



V. EXACT COVER 3 

A paradigmatic example for an NP-complete problem 
is exact cover 3: In this problem, one wants to find a 
bit-string Z\Z% . . . z n —\Z n , where the n bits z a € {0,1} 
must satisfy m constraints (clauses). For general 3- 
satisfiability (3-SAT), each of these clauses involves three 
bits a, /3, 7 £ {1, . . . ,n} and for the specific exact cover 
5-problcm every clause is defined by the constraint 



Z/3 



1 . 



(30) 



which has to be satisfied for every triple (a, j3, 7). A 
Hamiltonian encoding the solution to the exact cover 3- 
problem in its ground state can be defined by performing 
a sum over the (positive semidefinite) single-clause penal- 
ties (compare also [I3,|34j]) 



H F = 



c=l 



(*£ + • 



-1)* 



(31) 



Note that this problem Hamiltonian is slightly different 
from the original approach to exact cover 3 Q, which 
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assigns a fixed energy penalty to each violated clause 
and thus involves three-qubit interactions [26| . but has 
the same ground state (for satisfiable problems) - only 
some of the excitation energies differ. The above ansatz 
has the advantage that two-qubit operations suffice for its 
implementation, i.e., it is quadratic in the Pauli matrices. 
Inserting z,- = i [1 — of] one obtains 

H F = ml - + E ' (32) 

i i,j 

where > denotes the number of clauses involving the 
i th bit and riij > the number of clauses involving both 
bits i and j. 

For exact cover 3 our initial Hamiltonian (|29j) reads 

n 

B? = -Ef I*?*? + ■ ( 33 ) 

Note that due to the relations 2ni = n,j and 
3m = n ii the two-bit matrix of the final Hamilto- 
nian defines an exact cover 5-problem completely. 

In contrast, the original approach 5] employed the 
straight interpolation scheme with the initial Hamilto- 
nian 

^conventional = ^ [1 - of] , (34) 

i 

which has the unique ground state ([3]) and no broken 
symmetry. In addition, it only contains the single-bit 
structure of the final Hamiltonian. 

The Schrodinger equation is invariant under simulta- 
neous transformations of time and energy, such that we 
need to compare the energy scales of our modified ap- 
proach with the energy scales of the conventional ansatz 
[38| . The maximum energy spread of the quadratic final 
Hamiltonian (|32[) is AE max — 4m, which is only a fac- 
tor of four larger than the conventional final Hamiltonian 
Q . Likewise, the maximum energy of our initial Hamil- 
tonian (f3"3"]l can be upper bounded by Ai? max < 6m (the 
same bound would apply for the XYZ-network) , which is 
only a factor of two larger than E max = 3m for the con- 
ventional initial Hamiltonian (j3~4"|) . Therefore, the energy 
resources required by our modified algorithm are at most 
a factor of four larger than in the conventional approach. 

Regarding the unknown Hamming weight of the so- 
lution, we found for exact cover 3 that it is typically 
situated around A « n/3, which drastically reduces the 
polynomial overhead generated by trying every possible 
value of A, see Sec. HVl 

A. Hard exact cover 3 problems 

In the following, we will numerically compare the ef- 
fect of choosing the initial Hamiltonians (|3"3")l or (|34p on 
the performance of adiabatic quantum computation. Of 



course, this performance will depend on the selection of 
clauses in general. In order to select sufficiently hard 
problems, we only consider exact cover 3 problems with 
just an unique satisfying agreement, as was also done in 
the original study [5|. However, there is evidence that 
this restriction is not sufficient yet for ensuring the high- 
est computational complexity. For general 3-SAT, the 
number of instances with a unique satisfying agreement 
in combination with few clauses becomes exponentially 
rare among all instances with a unique satisfying agree- 
ment [24]. Also for exact cover 3 random problem in- 
stances with a unique satisfying agreement include many 
simple sub-problems, which have to be sorted out in or- 
der to find hard instances [35[ . Classically, the transition 
[39] ] from satisfiable to non-satisfiable problems for exact 
cover 3 suggests that problems with rather few clauses 
(m « 0.62n) have the highest computational complexity. 

In order to select problems near that classical phase 
transition in the small qubit range that is accessible to 
us, we have randomly generated two sets of 100 problem 
instances for each qubit number n: One problem set in- 
cluded only exact cover 3 instances that have a unique 
satisfying solution and an arbitrary number of clauses 
(as considered in 0]). Another problem class did include 
exact cover 3 instances with a unique solution and few 
clauses, specifically m < round(2n,/3). Note that this re- 
striction explains the slight triple clustering of datapoints 
around qubit numbers divisible by 3 in figures [H [1JJ [TTl 
[T2l and [121 Problem instances were generated similarly 
to Q by successively adding random clauses until there 
was only one solution left. In the first problem set (anal- 
ogous to [j| ) the result was discarded (and the procedure 
started over) whenever the problem became unsatisfi- 
able. In the second problem set we discarded the problem 
whenever either it became unsatisfiable or when the num- 
ber of clauses exceeded the boundary m < round(2n/3). 
Note that generating the latter problem set took a lot 
more effort than generating the first, which indicates that 
problems with a unique solution but few constraints are 
quite rare among all problems with a unique satisfying 
agreement - at least for the small qubit range accessible 
to us. 



B. Adiabatic Runtime 

Since analytic solutions are unfortunately not avail- 
able, we have to compare the performance of the con- 
ventional approach (|34|) with our modified proposal (|33|) 
numerically. To this end, we first compute the adiabatic 
runtime T that is necessary to obtain a fixed final fidelity 
of 1/8 for both schemes, see also [27j . In the conventional 
approach (|34[) . the quantum state was initialized with 
the state \S) in whereas it has been initialized with 
the normalized projection ([28)) of IS*) onto the correct 
Hamming subspace in our scheme. For a range up to 20 
qubits, the full Schrodinger equation had been integrated 
using a fourth order Runge-Kutta scheme [401 ] with an 
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FIG. 8: [Color Online] Runtime necessary to yield a final 
ground state fidelity between 0.12 and 0.13. In all cases the 
new algorithm (orange) performed in average better than the 
conventional approach (black). Hollow symbols correspond 
to especially hard instances of the exact cover 3 problem that 
do not only have a unique solution but also few clauses, in 
this case m < round(2n/3). (This leads to a triple cluster- 
ing of the hollow symbols around qubit numbers divisible by 
3.) Dashed lines show polynomial fits, whereas dotted lines 
display exponential fits (values given in brackets). 



adaptive step-size and varying run-times T, until an ac- 
ceptable success probability in the final state was found, 
see figure [SJ A considerably improved algorithmic per- 
formance was found - especially for the hard exact cover 
3 problems that have rather few clauses. For these hard 
problems, the conventional scheme performed in average 
significantly worse, whereas the new scheme did perform 
approximately similar on both problem classes. 

However, it should be emphasized that the perfor- 
mance on a specific problem may deviate significantly 
from the average behavior: For example, we have also 
found some problems where the conventional scheme per- 
formed better than the new algorithm. The median plot- 
ted in figure [8] for reasons of visibility of error bars and 
compatibility with || is hardly sensitive to such rare in- 
stances. The worst case runtime we find (not shown) is 
not even a monotonously growing function of the num- 
ber of qubits. For a given problem, it is therefore always 
practical to apply both quantum algorithms, since the 
solution can be tested in polynomial time. 



FIG. 9: [Color Online] Lower part of the spectrum for a typ- 
ical exact cover 3 instance with 10 qubits. After the position 
of the minimum gap was determined by a minimization algo- 
rithm, parabolae (dashed lines) were fitted to the lowest two 
eigenvalues. 
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FIG. 10: [Color Online] Value of the energies of ground state 
(solid symbols) and the first excited state (hollow) at the po- 
sition of the minimum gap for the conventional interpolation 
scheme (black) and the new approach (red) versus the num- 
ber of qubits. Each symbol corresponds to the median of 
100 random instances with a unique satisfying agreement and 
m < 2/3n clauses and error bars give the 99% confidence 
interval on the median. Dashed lines display linear fits (fit 
parameters given in the legend in brackets) to the ground 
state energies and demonstrate that the modified algorithm 
has only a moderate increase in the ground state energy. 



C. Spectral Properties 

In order to clarify the origin of the speed-up found 
for our scheme, we have also analyzed the lower part 
of the spectrum with the ARPACK package (37[: After 
determining the position of the critical point with a min- 
imization algorithm [47] applied to g{s) — E\(s) — Eq(s), 
a parabola was fitted to Eq(s) and to -Ei(s), see figure 
IH1 Since determining the lower part of the spectrum at a 



defined position is more efficient than integrating the full 
Schrodinger equation, we could extend our data range 
up to 23 qubits. In addition only the hard problem set 
with few clauses was considered here. The parameters of 
this parabola have been averaged over the 100 problem 
instances and are displayed in figures [TQl [Til an d Ell re- 
spectively. As one would expect, the scaling of the critical 
ground state energy is roughly linear for both algorithms, 
just the slope differs, see figure [TU] It is also visible from 
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FIG. 11: [Color Online] Value of the derivative of the lowest 
two energies at the position of the minimum gap. Invisibility 
of the hollow symbols demonstrates that an extremal point 
of the fundamental gap has been found with high accuracy. 
Color coding and statistics are analogous to figure [TOl 
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FIG. 12: [Color Online] Value of the curvature of the lowest 
two energies at the position of the minimum gap. Dashed 
lines represent polymial fits (parameters shown in brackets). 
Color coding and statistics are analogous to figure [TUl 



figure [TU] that the gap becomes smaller with increasing 
system size (as expected) and that the new approach has 
a larger minimum gap. The similar qualitative behavior 
holds also true for the average slope of the ground state 
energy at the critical point, see figure QT] However, for 
the curvature of the ground state energy - which is a 
direct marker for the order of the quantum phase tran- 
sition, compare subsection III Bl - we obtain significant 
differences between the two algorithms. Whereas for the 
conventional algorithm the magnitude of the curvature 
of the ground state (and the first excited state) at the 
critical point increases strongly, for the new algorithm 
this scaling is drastically reduced - though still existent. 
This provides some evidence that by choosing the initial 
Hamiltonian of the form (|33p with an approximately bro- 



ken symmetry and retaining the two-bit structure of the 
final Hamiltonian one may indeed improve the order of 
the phase transition for exact cover 3. 



D. Behavior of the fundamental gap 
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FIG. 13: [Color Online] Scaling of the value of the minimum 
gap versus the number of qubits. As one would expect, both 
scaling behaviors show a decrease with the system size. The 
vital question whether this scaling is polynomial (dashed fit 
lines) or exponential (dotted fit lines) cannot be answered 
with the limited numerical data. The minimum gap of our 
modified scheme is a factor of about 2.5 larger than in the 
conventional approach. Color coding and statistics are anal- 
ogous to figure [TOl 

Although the absolute value of the minimum gap is 
larger in our new algorithm than in the conventional ap- 
proach (see figure [T3")) , the ratio of the two values does 
not seem entirely sufficient for explaining the strong dif- 
ferences in the adiabatic runtime in figure [5J Recalling 
the results of section IIII1 we see that these results are 
compatible with the different behavior of the gap curva- 
ture, which can be deduced from figure [T2l and is shown 
explicitly in figurefTil This is evidence that the adiabatic 
runtime can be positively influenced not only by the min- 
imum fundamental gap alone but also by its curvature 
at the critical point. Note however that this is not too 
surprising, since for a smooth fundamental energy gap 
with a single minimum and that is bounded from below 
initially and finally, purely geometric arguments suggest 
that minimum gap and curvature at the minimum gap 
are related. 



E. Entropy of Entanglement 

Apart from the spectral properties of the Hamilto- 
nian, entanglement is another very useful concept for the 
understanding of quantum phase transitions as well as 
quantum computing 0, [4l| . If a quantum system is in 
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FIG. 14: [Color Online] Scaling of the fundamental gap cur- 
vature at the minimum gap versus the number of qubits. The 
curvature of our modified algorithm grows much slower as a 
function of n than with the conventional approach. Color 
coding and statistics are analogous to figure [TOl 



a pure state p = \^>) (>F|, the entanglement between two 
subsystems can be quantified by the von-Neumann en- 
tropy 

Se = S(pi) = -Tri {pilog 2 pi} 

= S(p 2 ) = -Tr 2 {p 2 log 2 p 2 } , (35) 

of the reduced density matrix of either subsystem 

Pi = Tr 2 {p} , p 2 = Tri {p} , (36) 

where Tri and Tr 2 denote the average over the degrees of 
freedom of subsystem 1 and 2, respectively. The quantity 
<Se is called entanglement entropy and has the property 
to vanish if and only if system and reservoir are not en- 
tangled, i.e., S E = <^ |*) = ® |* 2 >- It is also 
bounded by the smaller number of qubits in cither sub- 
system Se < min{ni,n 2 }. 

For general (sequential and adiabatic) quantum algo- 
rithms, it can be shown that an exponential speed-up 
is only possible if the entanglement Se contained in the 
system (at some point of the calculation) grows suffi- 
ciently fast with system size [42]. On the other hand, a 
nonanalytic (e.g., diverging) behaviour of the entangle- 
ment Se at the critical point (in the infinite size limit) 
is also a typical feature of quantum phase transitions. 
Let us illustrate these results by means of a few exam- 
ples: The Ising model discussed in subsection flTB] (which 
is analytically solvable) displays a logarithmic scaling of 
the entanglement entropy near the critical point: The 
ground state entanglement entropy of a block of L spins 
scales mildly as Se ~ \nL, see e.g., [43| and references 
therein. In contrast to spin systems in higher spatial 
dimensions with a typical power-law scaling Se ~ n p , 
for example, this logarithmic growth rate would not be 
sufficient for an exponential speed-up. Shor's algorithm 



contains order-finding [l|, where the entanglement en- 
tropy between source and target register scales approx- 
imately linearly with with the system size n [4lj ] . The 
same scaling has been found (numerically) for adiabatic 
algorithms for exact cover 3 by calculating the entangle- 
ment entropy for the first n/2 spins with the rest [411 ]. 
In comparison, the adiabatic Grover algorithm (with the 
same bipartition chosen) does also assume a maximum 
value of the entanglement entropy at the critical point - 
but this value remains constrained by one for all system 
sizes fill ]. 

Here, we have numerically analyzed the behavior of the 
entanglement entropy of the instantaneous ground state. 
Imagining the n qubits as being lined up in a chain, we 
calculated the entanglement entropy between the sub- 
systems formed by the first n\ qubits and the remaining 
n 2 = n — ri\ qubits. Numerically, the entanglement en- 
tropy was found by calculating the instantaneous ground 
state l^oW) of the system [37| and then determining the 
reduced density matrix p\. Due to S(pi) — S(p2) the 
largest reduced density matrices for n = 10 qubits have 
N 2 = 32 2 entries and can therefore be directly diagonal- 
ized [401 ] . The n — 1 values for the entanglement entropy 
resulting from the different partitions of an n-bit chain 
have simply been averaged. This resulting average was 
then again averaged over 100 (hard) problem instances 
of exact cover 3. 

As one might have expected, the curves for the conven- 
tional scheme and our proposal differ strongly, see figure 
1151 For the conventional scheme, the entanglement en- 
tropy vanishes at the beginning and at the end of the 
calculation and possesses a pronounced peak at the criti- 
cal point. In our algorithm, this peak is smeared out and 
we already start off with a relatively large entropy - the 
initial state in ([2"5|) is entangled, for a discussion of the 
required resources see subsection llV Al Moreover, the en- 
tropy is significantly larger in our algorithm throughout 
the interpolation. Both observations (i.e., the fact the 
the entropy is larger and varies slower) could be inter- 
preted as indications for the increased algorithmic per- 
formance - but one should bear in mind that the ex- 
act relation between the entanglement entropy and the 
achievable speed-up is not fully understood yet. 

Increasing the number of qubits, we do also see the 
approximately linear scaling of the maximum entangle- 
ment entropy (in the qubit range that is accessible to us) 
observed in the literature 13411. 



VI. NUMBER FACTORIZATION 

Apart from exact cover 3, there are many other prob- 
lems whose solution can be encoded in the ground state 
of a quadratic Hamiltonian in (|25[) . As one example, let 
us discuss the factoring problem (which is in NP, but 
not believed to be NP-complete) , a further example is 
outlined in the appendix IVIIII 

Given a large number lo, it is in general very difficult to 
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FIG. 15: [Color Online] Behavior of the average entangle- 
ment entropy for a system of 10 qubits. The shaded regions 
correspond to one standard deviation and the lines show the 
average over 100 problem instances. 



answer the question which numbers a and b fulfill lu = ab. 
The product of a /c-bit number a and another (n — A:)-bit 
number b has either (n — 1) or n bits. In the following, 
we assume the latter case (the other option can be easily 
adapted by using lu± = 0) and restrict our considerations 
to a Hilbert space with n qubits where the basis states 
can be written as |^) = \a%) . . . \a,k) |6i) ■ ■ • \b n -k)- Of 
course, even for problems with a unique solution (e.g., bi- 
primes with k ^ n/2), the partition (fc, n — k) will usually 
not be known in advance. However, the associated worst- 
case overhead of trying every possible value of k from 
n/2 to n grows only linearly as a function of n. A naive 
Hamiltonian encoding the solution as the ground state is 
then readily given by 



= (u>- ab) 2 



(37) 



where 



b = 



^a t 2 

1=1 

n—k 

J>2 



k-i 



1 



K 1 - "?) 



n-k-t 



1 = 1 



b e = -(l-o-l^) 



(38) 



Although it operates on a Hilbert space with dimension 
n and uses interactions between at most four qubits, the 
above Hamiltonian has the disadvantage that its cou- 
plings (operator pre-factors) cover an exponential range 
and thus its spectral range (ratio of largest and smallest 
eigenvalue) increases exponentially with the number of 
qubits n. In an experiment, an exponential fine-tuning 
of the couplings will probably be infeasible for the in- 
teresting case of large n. In addition, an exponentially 
increasing spectral width will be hard to realize without 
using exponential resources. 

This scaling of the spectral width could be improved by 



using the Hamming distance as a penalty in the Hamil- 
tonian 



1=1 



[u t + (-i) ut (aS) 



(39) 



where u>i denotes the l-th. bit of u, since here, the spec- 
tral width grows only linearly with the number of qubits. 
However, the computation of the £-th bit of the product 
ab requires simultaneous interactions of multiple qubits. 

In the following, we construct a Hamiltonian that 
has at most two-qubit interactions and simultaneously 
a spectral width that scales only polynomially with the 
number of qubits. In addition, we demand that the local 
coupling constants should cover a finite range only. 



A. Factorization Equations 

To introduce our approach, let us consider as a generic 
example the case n = 10 and k — 6. Then, the usual 
method of multiplying two binary numbers can be de- 
composed into n equations according to Table Q] such as 



ai&4 + a 2 b 3 + a 3 b 2 + a 4 6i + z 65 + z 75 
- (u> 5 + 2z 54 + 4z 53 ) 



, (40) 



where the binary variables represent the carries from 
columns i towards j (the above equation corresponds to 
the fifth column). The complete set of equations like 
(|40| are also called factorization equations, where a linear 
(plus logarithmic corrections) amount of carry variables 
Zij is needed, see 0| for a detailed discussion. In an 
adiabatic quantum computer, these carry variables would 
have to be represented by additional qubits. 
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TABLE I: Conventional method for multiplying the k = 6 bit 
number a and the n — k — 4 bit number 6. The result uj = ab 
is obtained by summing the partial products aibj and carry 
bits (from column i to column j), which gives rise to n 
factorization equations. Note that for this example, at most 
two carry bits are necessary in each column. 

Obviously, factoring the number lu corresponds 
to finding a solution for the binary variables 



(ai, . . 
tions. 



, a kl bi, 



,b n ^k,Zij) to the factoring equa- 



A simple idea is to assign a positive penalty to 
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each violated factoring equation via 



E^ 2 



(41) 



where Ej = expresses the ith factoring equation in nor- 
mal form (with zero on the right hand side). Note how- 
ever, that by naively using the square of the factoriza- 
tion equations we have generated four-qubit-interactions 
in the Hamiltonian. The maximum penalty originating 
from a single equation scales as 0{(n — k) 2 }, such that 
the complete spectral width of this Hamiltonian scales 
in the worst case as 0{n 3 } (and so does the number of 
required four-qubit interactions). Although only a lin- 
ear number of auxiliary qubits is needed, it may be ex- 
perimentally difficult to realize so many different four- 
qubit interactions. In addition, our approach to generate 
higher-order phase transitions cannot be directly applied, 
since the Hamiltonian does not have the standard form 
of (J25l). 

In order to obtain only two-qubit interactions, further 
modifications are necessary: Consider the case where 
only a single product of binary variables A and B is in- 
volved in each factoring equation 



E = AB + S = , 



(42) 



where A and B are single bit variables and S = ■ OjCj 
is a sum of single bit variables with ay G Z. Then, 
the naive penalty for violating the equation (|4"2")) . which 
would be given by Pi = [AB + S] , can be replaced by 



Pe = 2 



S 



(43) 



in the sense that both penalties vanish if and 
only if equation (|4"2"| is obeyed [which implies 
(AB = A S = 0) V (AB = 1 A S = -1)] and are larger 
or equal to one otherwise. This modified penalty has the 
advantage that it involves only quadratic interactions be- 
tween the qubits. 

The task one is left with is to rewrite the multiplication 
table with the help of further ancilla variables such that 
only factorization equations of the type (|4"2")) occur. For 
the generic example, this is demonstrated in table ITT1 

The general structure of the equations is for 1 < i < k 
and 1 < j < n — k is given by 



a t bj 



Sij T Zi 



Si-i 



1,3-1 



2z 



(44) 



where at the boundaries (marked by invalid indices) the 
equations are closed by 



z 0,j 




Sifi 


= Wt, 




= ^k+j , 


Si ,n — k 


z i,n—k 


Z k ,j 


= o, 


Sln-k-1 


= 0. 



(45) 
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TABLE II: Multiplication table for multiplying the 6 bit num- 
ber a and the 4 bit number b using only two-qubit interac- 
tions. In each row, the first line contains the partial product 
sums, the second line the partial products, and the last line 
the carry variables. Note that the meaning of the indices 
for the auxiliary variables has changed in comparison to ta- 
ble U Apart from the boundaries, the binary variables have 



tO fulfill Qibj + Si 



= Si 



+1,3 



-1 + 2zi-ij (see the text for 



explanations). For each of the (n — k — 1 = 3) horizontal lines 
that have been inserted, k = 6 partial product variables Sij 
and (k — 1 = 5) carry variables 2y have to be included. The 
first partial product variable (second row) always vanishes. 
Note also that the first partial product variable in each row 
is equal to the carry variable in the row above, which has 
therefore been removed. 



With the suggested replacement in (|43|) , this leads to the 
total penalty Hamiltonian 



(3) _ 



k 71 — k 

EE] 2 

»=i j=i 



- (&i +bj - -) + S i3 -r -, , 



Si+XJ— 1 — 2£;_ij 



(46) 



where all operators act on distinct single qubits (for 
example on a two-dimensional lattice arrangement as 
Sij = (1 — <7?-)/2). The above Hamiltonian has the 
ground state (an appropriate qubit ordering assumed) 



Ol 



.a fc )|6i...& n _ fe )|{^-})|{ % }) , (47) 



where to = ab is the sought-after factorization. The num- 
ber of necessary auxiliary variables can be calculated as 
follows: In order to multiply a fc-bit number a and a 
(n — fc)-bit number 6, at most k(n — k) partial prod- 
ucts are required. For separating these into sums with a 
single product of two binary variables, (n — k — 1) hori- 
zontal lines have to be inserted, compare table OH Each 
of these lines requires k partial sum variables Sij (except 
in the second row, where the first Si,n-fc-i vanishes) and 
(k — 1) carry variables Zij . Note that the bottom (k — 1) 
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carry variables can be associated with the vanishing top 
(k — 1) carry variables. Therefore, the total number of 
auxiliary variables is given by (2k — l)(n — k — 1) — 1. 
Together with the n bits required for the factors a and 
b, at most n — 1 + (2k — l)(n — k — 1) qubits are re- 
quired for finding the k and n — k bit factors of an n-bit 
number. On an adiabatic quantum computer, this im- 
plies a quadratic overhead in the number of qubits. Since 
the number of variables in each factoring equation is al- 
ways smaller or equal to 6, the total number of quadratic 
interactions scales with the number of equations, i. e., 
quadratically in n. For each equation in (|46p . the largest 
possible penalty is 21, which leads to a quadratical scal- 
ing of the spectral width. Note also that the necessary 
coupling strength between different qubits ranges (inde- 
pendent of the problem size) from 1 to 8. 



B. Numerical Study 

Similar to satisfiability problems - where those with 
a unique solution are believed to belong to the hardest 
problems [f| (both classically and in adiabatic quantum 
algorithms), bi-primes (products of two prime numbers) 
possess a unique factorization and the classical hardness 
of bi-prime factorization is used in many cryptography 
protocols. Except in the case k = n/2 (where the solution 
becomes two-fold degenerate due to u> — ab — ba), the 
Hamiltonians have a unique ground state. In our case, 
biprime factorization constitutes a well-defined (and also 
classically relevant) problem class, for which the algo- 
rithmic performance of an adiabatic quantum algorithm 
using (|46[) as the problem Hamiltonian would be inter- 
esting. 

In order to reduce the problem complexity of (|46[) , one 
can use the fact that both prime factors are odd in non- 
trivial cases and that there exists a minimum size of the 
prime factors to generate an n-bit number, i.e., that both 
first and last bits of the prime factors are set to 1, see ta- 
ble [ITT] for the generic example. Under these conditions, 
the first row generates k — 1 nontrivial equations that in- 
volve no product between different qubits at all. In these, 
the Sij variables can be eliminated without changing the 
structure of the equations. Together with the (n — 4) 
bits required to store the unknown bits of the prime fac- 
tors, one arrives at [2k(n — k — 1) — 3] variables to find 
the odd (k, n — k) factors of the n-bit number u>. Finally 
the two-qubit Hamiltonian can be cast into the form 

"tot "tot "tot 

H F = h+j2 h i a ! + 2 J2 J2 h ij°ivh ( 48 ) 

i— 1 i— 1 j=i J r\ 

where ntot denotes the total number of qubits and 
h,hi,hij £ M with hij — hji and ha = 0. It is easy 
to see that the coefficients h, hi , hij can be conveniently 
extracted from eqn. PB) by computing traces of , 
i4 3) , and erf a j Hp \ respectively. We have analyzed 



the algorithmic performance for a linear interpolation for 
three different initial Hamiltonians 

"tot -, 

i=l 

"tot I , I 

H7 = £ ™ [2 - « - crV a]] , 

"tot I , I 

H ,y Z = £ KM [x _ gXgX _ a y a y _ ^ (4Q) 

Note that the first Hamiltonian has the ground state \S) 
in ([3]) and no apparent symmetry of the final Hamilto- 
nian is respected. The latter two Hamiltonians commute 
with the Hamming-weight operator (|26[) and the same 
discussion as in section [V] applies. However, there are 
also some crucial differences: 

• In contrast to exact cover 3, the numbers hij in the 
quadratic decomposition (|48p are not necessarily 
positive. 

• In addition, unlike exact cover 3, the two-bit inter- 
actions hij do not completely determine the prob- 
lem, such that the Hamiltonians in (|4"5|) can prob- 
ably be further improved. 
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TABLE III: Multiplication table for multiplying the odd 6 
bit number a and the odd 4 bit number b. The (k — 1) 
Sij -variables in the top row of table [Hi have been eliminated 
without changing the structure of the equations. One equa- 
tion (top right) becomes trivial, whereas (from top to bottom) 
fc — 1 + 2 (n — k — 2) + k = 2n~ 5 equations emerge that involve 
no products between different qubits at all. 

Since we can with moderate effort simulate the evo- 
lution of systems with 0(20) qubits, we can numerically 
access the (n, n— k) factoring partitions displayed in table 
IIVI For the bi-primes within this range, we have deter- 
mined the minimum fundamental gap during the linear 
interpolation (3?| , see figure [16j It is visible that again 
the H^ y Hamiltonian is superior to the other choices in 
(|4^|) - at least for the small sample problems under con- 
sideration. 
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Interestingly, the third Hamiltonian in (|49[) leads in 
some cases to an even smaller minimum gap than with 
the conventional choice. We have also observed this for 
rare instances of EC3. As discussed in subsection IIVB1 
we conjecture additional nearly conserved quantities to 
hamper the adiabatic evolution towards the solution in 
case of the H^ yi initial Hamiltonian. 



n to t 


bi-primes 


partitions (n, n — k) 


5 


33;39 


(4,2);(4,2) 


7 


51;57;69;87;93 


(5,2);(5,2);(5,2);(5,2);(5,2) 


9 


25;35;49;111;123 
129;141;159;177;183 


(3,3);(3,3);(3,3);(6,2);(6,2) 
(6,2);(6,2);(6,2);(6,2);(6,2) 


13 


55;65;77;91 


(4,3);(4,3);(4,3);(4,3) 


17 


85;95;115;119;133 
145;155;161;203;217 


(5,3);(5,3);(5,3);(5,3);(5,3) 
(5,3);(5,3);(5,3);(5,3);(5,3) 


21 


121;143;169 


(4,4);(4,4);(4,4) 



TABLE IV: Bi-primes and partitions accessible with different 
numbers of qubits (leftmost column). 



O 
& 



BO 

■B 

T3 



• 2 qubil final, X initial 
□ 2 qubit final, XY initial 
A 2 qubit final, XYZ initial 



4 6 8 10 12 14 16 18 20 22 

total number of qubits 

FIG. 16: [Color Online] Minimum gap during adiabatic inter- 
polation for number factorization versus the total number of 
qubits used. For problems with a doubly degenerate ground 
state, the minimum gap between the lowest and the second 
excited state has been calculated. Each data set has the same 
number of data points per column (some symbols lie on top 
of each other). 



(cf. Fig. [4]), we conjecture that they are advantageous 
with respect to adiabatic quantum computation in the 
sense that they allow much shorter run-times. 

Based on this physical intuition, we designed a mod- 
ified adiabatic quantum algorithm for the NP-completc 
problem exact cover 3 in analogy to a symmetry-restoring 
quantum phase transition and found numerically that it 
indeed yields an improved performance in comparison 
with the conventional scheme. Even though the infinite- 
size limit of these adiabatic quantum algorithm may not 
be as well-defined as that of the Ising model (probably 
they should be classified somewhere between first and 
second order), the spectral characteristics indicate that 
our alternative adiabatic quantum algorithm has a higher 
order than the conventional scheme. 

As a closely related point, we observed that the adi- 
abatic run-time (of constant-speed interpolations) is not 
just determined by the value of the minimum gap but 
also by the curvature of the energy levels at the critical 
point. In fact, the better average scaling behavior of our 
alternative algorithm can probably be attributed mostly 
to the fact that the energy levels are far less curved than 
in the conventional scheme. 

Apart from the spectral characteristics, entanglement 
is another useful concept for quantum algorithms as well 
as quantum phase transitions. We found that the ground- 
state entanglement entropy in our modified adiabatic 
quantum algorithm is much larger than that in the con- 
ventional scheme, which might be connected with the 
advantages of higher-order phase transitions. However, 
these advantages may also go along with some drawbacks: 
for the prototypical Grover and Ising models, we found 
that second-order phase transitions seem to be more vul- 
nerable to decoherence than those of first order [13, EH • 

Finally, even though the run-time scaling (polynomial 
versus exponential) is not clear, the results obtained so 
far strongly suggest that adiabatic quantum algorithm 
algorithm can solve NP problems much faster than the 
Grover search routine (with a quadratic speed-up). The 
methods proposed in this article can easily be applied to 
other computationally relevant problems which can be 
encoded into the ground states of Hamiltonians that are 
quadratic in the Pauli matrices - such as factoring (see 
also the Appendix IVIII[) . 
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VIII. APPENDIX: ENCRYPTION WITH TEA 

Although factoring numbers is also used in encryption 
algorithms, it is by no means the only way to encode 
or decode information. One further example is the Tiny 
Encryption Algorithm (TEA). The C-code for TEA [Iff 
is publically available and the algorithm has not been 
broken - i.e., there is no known classical algorithm which 
finds the key significantly faster (with input and output 
given) than the simple brute-force search. In its standard 
form, it encodes 64 input bits using a 128 bit key. En- 
coding is done by an iterative procedure that performs 
- depending on the values of the key bits - different op- 
erations on the input bits. With knowledge of the key, 
these operations are easily invertible, whereas without 
the key the inversion of the encryption becomes expo- 
nentially complex. Therefore, extracting the key from a 
given input and its encrypted output would open a way 
to decode all information that has been encrypted with 
this key. 

Here we will demonstrate that in principle it is possi- 
ble - if input and output are known (e.g., for a specific 
communication) - to encode the search for the key in the 
ground state of a quadratic Hamiltonian of the form (f25|) . 
With the input given by the n — 32 bit numbers y 1 and 
z 1 and the output by the 32-bit numbers y l and z , the 
TEA algorithm establishes a relation between these and 
the encryption key numbers k a , k , k c , kd by forming the 
sequence 

y l+1 = [S^(z l ) + k a ] © [z l + «*] © [S+(z l ) + k b ] , 
z t+1 = [Sztf) + ko] © [y l + «'] © [S£ (if) + k d ] , 



where the index i runs from 1 to £ — 1. Here © de- 
notes bitwise addition (XOR), defines left- and right- 
bit-shifting by a digits, and s l are some fixed (known) 
numbers given by the algorithm. In its standard form, 
the algorithm uses i — 32 iterations. Obviously, in each 
of the three different operations used - adding of two 
numbers, bit shifting and bitwise addition (XOR) - al- 
ways one of the two operands is defined by the key or by 
the encryption algorithm. Therefore, knowledge of input, 
output and the key enables one to reverse the process (of 
course, with knowledge of the encryption algorithm as- 
sumed). From the previous section it is evident that the 
addition of two n-bit numbers [none of which is known - 
e.g., S^(z l ) + k a ] can be encoded in the ground state of 
a quadratic Hamiltonian using just n ancilla bits. Using 
the replacement (|4"3"]) . the same is possible for the XOR 
operation. Note that the bit shifting used does not pose 
any (theoretical) problem, since it just changes the in- 
dex of qubits which have to interact. In each iteration 
of (|50[) one has 4 XOR operations and 4 non-trivial ad- 
ditions (where both addends are unknown). Therefore, 
a quadratic Hamiltonian would not only operate on the 
key bits k a ,kb,k c , kd but also on the intermediate results 
y l ,z % as well as the ancilla qubits. In each iteration, 
0{n} auxiliary qubits are required, such that an overall 
number of 0{n£} bits will suffice to encode the search 
for the key in the ground state of a quadratic Hamilto- 
nian. Ergo, similar to factoring, we obtain a quadratic 
overhead and thus one might expect that an adiabatic 
quantum algorithm might find the key much faster than 
the Grover (brute-force) search. 
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the binary numbers by adding terms of the form z 2 — zt. 
Although the minimum of a positive definite quadratic 
form in n dimensions can easily be found in the set of 
real numbers with a standard minimization algorithm, 
is must be kept in mind that the above-posed problem 
implies the side constraints that the numbers have to be 
binaries. These side constraints have to be expressed by 
higher-order polynomials [such as e.g., J^a^l — Xi) 2 ], 
which lead to the existence of many local minima and 
classical minimization algorithms will have difficulties to 
find the global minimum. 

[47] Of course, for the determination of the eigen-energies of 
the new approach we did only consider the relevant Ham- 
ming subspace. 



